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CHAPTER 1 
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In a two-dimensional version of the modified Hasegawa-Wakatani (HW) 
model, which describes electrostatic resistive drift wave turbulence, the 
resistive coupling between vorticity and density does not act on the 
zonal components {ky = 0). It is therefore necessary to modify the HW 
model to treat the zonal components properly. The modified equations 
are solved numerically, and visualization and analysis of the solutions 
show generation of stable zonal fiows, through conversion of turbulent 
kinetic energy, and the consequent turbulence and transport suppres- 
sion. It is demonstrated by comparison that the modification is essential 
for generation of zonal flows. 

1. Introduction 

In quasi two-dimensional (2D) plasma and fluid flows the energy flux from 
small scale turbulent modes toward lower wavenumber modes can domi- 
nate the classical Kolmogorov cascade to dissipative scales, with the result 
that energy can accumulate in large scale coherent structures. Zonal flows 
in planetary atmospheres and in magnetically conflned fusion plasmas are 
well-known examples of such coherent structures. Quasi two-dimensional 
fluid systems in which turbulent activities and coherent structures interact 
can undergo a spontaneous transition to a turbulence-suppressed regime. 
In plasmas such transitions dramatically enhance the conflnement and are 
known as L-H or confinement transitions. From theoretical and experimen- 
tal works the importance of shear or zonal flows for suppression of cross-fleld 
transport and conflnement improvement is now widely appreciated. 

Several low-dimensional dynamical models, comprised of a small num- 
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ber of coupled ordinary differential equations, have been proposed to de- 
scribe and predict the L~H transition '-' . Ball et al. have analyzed a 
three- variable model using bifurcation and singularity theoriei^. The model 
is based on the reduced resistive magnetohydrodynamic equations with 
the electrostatic approximation, and describes the pressure-gradient-driven 
turbulence-shear flow energetics. This approach using low-dimensional 
modeling greatly simplifies the problem, and when validated against simu- 
lated or real experimental data, will provide an economical tool to predict 
transitions over the parameter space. 

In this work we report the results of numerical simulations that both 
complement the low-dimensional modeling results and raise some interest- 
ing issues in their own right. We focus on a model for electrostatic resistive 
drift wave turbulence, the Hasegawa-Wakatani (HW) modePl, and solve 
the equations by direct numerical simulation in 2D slab geometry. The 
HW model has been widely used to investigate anomalous edge transport 
due to coUisional drift waves^. Moreover, self-organization of a shear flow 
has been shown by numerical simulation of the HW model in cylindrical 
geometrj0. Thus we consider the HW model is a good starting point for 
studying self-consistent turbulence-shear flow interactions, even though it 
does not describe physics that can be important in specific situations, such 
as magnetic curvature, magnetic shear, and electromagnetic effect. 

2. Modified Hasegawa-Wakatani Model 

The physical setting of the HW model may be considered as the edge region 
of a tokamak plasma of nonuniform density hq = no(a:) and in a constant 
equilibrium magnetic field B = Bq\7z. Following the drift wave orderinj^, 
the density n = jiq + ni and the electrostatic potential ip perpendicular 
to the magnetic field are governed by the continuity equation for ions or 
electrons and the ion vorticity equation, 

^d?^^^-^°a;^- 

where ~ (d/dx, d/dy)^ , d/dt ~ d/dt + V e ■ is the E x B convec- 
tive derivative (Ve = ^ Vj^^s x Wz/Bq, E = —V±(p), ni is the ion mass, 
jz is the current density in the direction of the magnetic field. The conti- 
nuity equation ^ can refer to ions and electrons because V • j = under 
the quasineutral condition, and ([2]) holds because the current density is 
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divergence-free. Since the ion inertia is negligible in the parallel direction 
(z), the parallel current is determined by the Ohm's law, 

E+—Vp,^Tjj. (3) 

en 

If the parallel heat conductivity is sufficiently large, the electrons may be 
treated as isothermal: Pc = nT^ (p is the pressure, T is the temperature, 
and subscript e refers to electrons.) This gives the parallel current as 

3z = — 1'^" ■ 

7] oz \ e J 

If we eliminate from ([1]), Q and normalize variables as 

x/pb-^x, ujcit t, eip/Tc-^ifi, ni/nn—^n, (5) 



-1 



where Wci = eBo/m is the ion cyclotron frequency, and ps = ^Tc/mui^ 
is the ion sound Larmor radius, we finally obtain the resistive drift wave 
equations known as the Hasegawa-Wakatani (HW) modePl, 

^C + {v',C} = «(^-^)-^cV'C (6) 

—n + {ip,n}^a{ip-n)-K—-DnV'^n, (7) 

where {a,b} = {da/dx){db/dy) — {da/dy){dh/dx) is the Poisson bracket, 
= d"^ jdx^ + d"^ jdy^ is the 2D Laplacian, = V^i^j is the vortic- 
ity. We omit _L, and use V for the 2D derivative. The dissipative terms 
with constant coefficients I?^ and I?„ have been included as adjuncts 
without derivation, for numerical stability. The background density is as- 
sumed to have an unchanging exponential profile: k, = —{d/dx)l'nno. 
a = — To/(r/noWcie^)9^/9z^ is the adiabaticity operator describing the par- 
allel electron response. In a 2D setting the coupling term operator a be- 
comes a constant coefficient, or parameter, by the replacement d/dz —>■ ifc^. 
This resistive coupling term must be treated carefully in a 2D model be- 
cause zonal components of fluctuations (the ky ^ kz ~ modes) do not 
contribute to the parallel currenllSi. Recalling that the tokamak edge turbu- 
lence is considered here, ky = should always coincide with kz = because 
any potential fluctuation on the flux surface is neutralized by parallel elec- 
tron motion. Let us define zonal and non-zonal components of a variable / 
as 

zonal: if) = J ^^y^ ^on-zonal: f = f - (/), (8) 
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where Ly is the periodic length in y, and remove the contribution by the 
zonal components in the resistive coupling term in ^ and ([7|). By sub- 
tracting the zonal components from the resistive coupling term a{ip — n) — > 
a{(p — n), we end up with the modified HW (MHW) equations, 

^^C + {^,C} = ai^-n)-DcV% (9) 

—n + {(f, n} = a{(p - n) - n- D„V^n. (10) 

ot ay 

The evolution of the zonal components can be extracted from © and PT]|) 
by averaging in the y direction: 

|(/) + |:(>.) = /^v^(/), (11) 

where / stands for C, and n, and D stands for the corresponding dissipation 
coefficients. 

The HW model spans two limits with respect to the adiabaticity 
parameter. In the adiabatic limit a ^ oo (coUisionless plasma), the 
non-zonal component of electron density obeys the Boltzmann relation 
fi = no(2) exp((^), and the equations are reduced to the Hasegawa-Mima 
equation^. In the hydrodynamic limit a ^ and the equations are decou- 
pled. Vorticity is determined by the 2D Navier-Stokes (NS) equation, and 
the density fluctuation is passively advected by the flow obtained from the 
NS equation. 

In the ideal limit (a = oo, Dc^ — Z?„ — 0) the modified HW system has 
two dynamical invariants, the energy E and the potential enstrophy W , 

E=^Jin' + \\/ip\')dx, W^^Jin~C?dx, (12) 

where da; — dxdy, which constrain the fluid motion. According to Kraich- 
nan's theory of 2D turbulencJ^, the net flux of enstrophy is downscale 
while that of energy is upscale. This inverse energy cascade is behind the 
development of large scale, stable coherent structures in a HW flow. 
Conservation laws are given by 

— r- = r„ — Fq — De, —— = r„ — Dw- (13) 
at at 
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(14) 



(15) 



(16) 



Dw= [DniV^nf + D^{V\f - {D„ + Dc)\7^nV\]dx. (17) 



These quantities constitute sources and sinks. As will be seen in the simu- 
lation results, they are mostly positive (r^ and De are positive definite), 
thus only r„ can act as a source. The energy absorbed from the background 
supplies the turbulent fluctuations through the drift wave instability. 

Note that the same conservation laws hold for the unmodified original 
HW (OHW) model except that Ta is defined by both zonal and non-zonal 
components; r^^^ = a J{n — ip)^dx. In the OHW model, the zonal modes 
as well as the non-zonal modes suffer the resistive dissipation. 

2.1. Linear Stability Analysis 

Since the zonal modes have linearly decaying solutions, we only consider 
the form gi(fcxa;+/cai/-wt) (-^^ ^ Qy Linearization of the equations ([9]) and 
([T0|) yields the dispersion relation, 

u;'^ + iuj(b+{l + p-^)k^D^)-ibuj^-ak^{P + p-^)Di;-k^p-^Dl = 0, (18) 

where we defined — k"^ + k^, b = a{l + k"^) / k'^ , the drift frequency 
(jj* = kyK/{l + k"^), and the Prandtl number Pr = D(^/Dn- Solutions to the 
dispersion relation are given by 



cr = A,ak'^{P + p-^)D^+^k^p-^Dl-{b+{l+p-^)k^Di^f , ta.ne^ 'ibujja. 



In the limit where Dq — Dn = 0, it is readily proved that one of the growth 
rate 7 = 3(a') is positive if bcu^, is finite, thus unstable. However, there exists 
a range of where the drift wave instability is suppressed. The stability 
threshold is given by 



5R(w) = ±-(t72 + 166V') ^ cos-, 
3(c^) = -^ b+{l + p-^)k^DcT{<J^ 



+ 



(19) 

166V) ^ sin ^ , (20) 




(21) 
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and is depicted in Fig. [T] The left panel shows the stability boundary in 
— K plane. If we enhance the drive by increasing k, the system becomes 
unstable. However, the instability is stabilized by increasing the dissipation. 
The stability threshold in — ky plane is shown in the right panel. We 
see that in a highly driven-dissipative system only low wavenumber modes 
are unstable. The stability boundary in parameter space is a region where 
interesting dynamics are expected to occur, such as bifurcations or sudden 
changes to a suppressed (or enhanced) turbulence regime. 

Figure [2] shows the dispersion relation for cases where — £)„ — 0. To 
provide a test of the simulation code, we plot growth rates obtained from 
numerical simulations together with the analytic curves. We can see that 
the growth rates obtained numerically agree very well with that calculated 
analytically. We also note that, in the parameter range plotted in Fig. [2] 
{a = I, K = 1), the most unstable mode is k^ ~ 0, /cy ~ 1. 

3. Simulation Results 

The HW equations are solved in a double periodic slab domain with box size 
(2L)^ = (27r/AA:)^ where the lowest wavenumber Afc = 0.15. The equations 
are discretized on 256 x 256 grid points by the finite difference inethod. 
Arakawa's method is used for evaluation of the Poisson bracketfl^. Time 
stepping algorithm is the third order explicit linear multistep method^. 

Since we are focusing in this work on how the modification ([9]), (|10|) 
influences nonlinearly saturated states, we fix the parameters to k = 1, 
Z?^ = 10^^, a — 1, and P,. = 1, and compare the results obtained using 
the MHW model with those computed from the OHW model. For these 
parameters the system is unstable for most wavenumbers. During a typical 
evolution, initial small amplitude perturbations grow linearly until the non- 
linear terms begin to dominate. Then the system arrives at a nonlinearly 
saturated state where the energy input r„ and output due to the resistivity 
Ta and the dissipations De.w balance. 

In Fig. \3\ we contrast the zonally elongated structure of the saturated 
electrostatic potential computed from the MHW model with the strong 
isotropic vortices in that from the OHW model. Time evolution of the 
kinetic energy — 1/2 J \ W(p\'^dx, and its partition to the zonal and the 
non-zonal components are shown in Fig. [H The saturated kinetic energy 
is not affected by the modification {E^ ^ 1 for both cases). In the OHW 
model, the zonal flow grows in the linear phase, as well as the other modes, 
up to a few percent of the kinetic energy, and saturates. On the other hand. 
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Fig. 1. Stability diagram of the MHW model. Left panel shows the stability thresholds 
in — K plane. The drift wave instability can be stabilized by strong dissipation. In 



the right panel, stability thresholds are plotted in Kx 
only some low wavenumber modes are unstable. 



- ky plane. For certain parameters, 



in the MHW model the zonal kinetic energy continues to grow after the 
linear phase, and dominates the kinetic energy. The kinetic energy contained 
in other modes decreases to a few percent of the total kinetic energy. In the 
original 2D HW model, the resistive coupling term is retained for the zonal 
modes, the effect of which is to prevent development of zonal flows. But 
since the zonal modes do not carry parallel currents it is clearly unphysical 
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Fig. 2. Dispersion relation of the dissipationless MHW model, a = 1, k = 1. 

to retain resistive action on them. Subtraction of the zonal components 
from the resistive couphng term is necessary to permit the generation of 
zonal flows. 




Fig. 3. Contour plots of saturated electrostatic potentials for the modified and the 
original HW models. Zonally elongated structure is clearly visible for MHW case. 

The density flux in x direction r„ (transport across the magnetic field) , 
together with the energy partition to the kinetic energy and the po- 
tential energy = 1/2 Jn^dx, is plotted in Fig. \5\ We observe that 
once the zonal flow is generated in the MHW model, the transport level is 
significantly suppressed. The transport suppression is mostly because the 
saturated potential energy (or amplitude of saturated density fluctuation) 
is reduced. The potential energy and the turbulence kinetic energy are con- 
verted into the zonal kinetic energy. By contrast the energy of the OHW 
model is almost equi-partitioned between the kinetic and potential energy. 
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Fig. 4. Time evolution of the kinetic energy, and its partition to the zonal and the non- 
zonal components. In the modified HW model, the zonal mode contains most of kinetic 
energy, while non-zonal turbulence contains most of the kinetic energy in the original 
HW model. 



The kinetic energy spectra averaged over the x or y direction for the 
MHW and the OHW models are shown in Fig. [6l The x [y) averaged kinetic 
energy spectra ^(j,) ) ^^re defined from the Fourier amplitude of the kinetic 
energy £^ hy 

^C"-) = 1^ / " £'^ik,,ky)dky, (22) 

J^y Jo 

1 r^- 

£fiky) = — £'^ik,,ky)dk,, (23) 

J<x Jo 

where Kx,Ky are the highest wavenumbers. The spectra of the modified 
model again show strong anisotropic structure whereas there is no marked 
difference in the original HW model. In the modified model, potential energy 
stored in the background density is converted into turbulent kinetic energy 
through the drift wave instability at ky ^ 1, kx = and then is distributed 
to smaller wavenumbers. The drift wave structure, which is elongated in 
the X direction, is break up into rather isotropic vortices after the nonlinear 
effect sets in, and those isotropic vortices merge in the y direction to produce 
the zonal flow. We can recognize this non-negligible inverse energy cascade 
in the y direction from a slight negative slope of £x{ky) spectrum in ky 
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Fig. 5. Time evolutions of the radial density transport and the kinetic and the potential 
energies for the modified and the original HW models. Once zonal flow is generated in 
MHW model, the turbulent fluctuation level and transport are significantly reduced. 



region. The y averaged spectrum £x{ky) shows the strong peak at the zonal 
wave number ~ 0.45. 



4. Conclusion 

We have performed nonlinear simulations of the 2D HW model. As sug- 
gested recentlj^^, the electron response parallel to the background magnetic 
field must be treated carefully in the 2D model. The model should be mod- 
ified to exclude the zonal {ky = 0) contribution from the resistive coupling 
term. By comparing the numerical results of the modified and the unmod- 
ified original HW models, we have revealed that a remarkable zonal flow 
structure in the nonlinearly saturated state is only observed in the modified 
model. Thus, the modification is crucial to the generation of the zonal fiow 
in this model. Time evolutions of the macroscopic quantities, such as the 
energies and fiuxes show that, after the zonal flow is built up by turbulent 
interaction, the generated zonal flow signiflcantly suppresses the turbulent 
fluctuation level and the cross-fleld density transport. 

The build up of the zonal flow and resulting transport suppression in- 
dicate bifurcation structure of the system. If we increase a parameter (say. 
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Fig. 6. The x and y averaged kinetic energy spectra for the MHW and the OHW 
models. The top two lines (solid line for £y{kx) and broken line for £x(ky)) for the OHW 
model are almost overlapped indicating isotropy. The middle two lines (dot-dashed line 
for £y(kx) and dotted line for Ex{ky)) for MHW show highly anisotropic structure in low 
k region. The energy injected at {k^, ky) = (0, 1) cascades inversely to the zonal mode 
of the wave number (0.45, 0). The bottom two series of symbols show the linear growth 
rates of modes for reference. 



strength of the hnear drive term k), the system may undergo sudden tran- 
sition from a high transport to a low transport regime. The state shown 
in this paper can be a bifurcated state. A systematic parameter study and 
comparison with the low-dimensional dynamical model are possible next 
steps. 
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